{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import preprocessing as pre\n",
    "import numpy as np\n",
    "from numpy import linalg as la\n",
    "import random\n",
    "import math\n",
    "import sklearn\n",
    "import scipy.io as sio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#X1(CENTRIST FEATURE):\n",
    "X11=np.load(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\dollar_bill\\CENTRIST_feature.npy\")\n",
    "X12=np.load(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\Faces\\CENTRIST_feature.npy\")\n",
    "X13=np.load(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\garfield\\CENTRIST_feature.npy\")\n",
    "X14=np.load(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\Motorbikes\\CENTRIST_feature.npy\")\n",
    "X15=np.load(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\snoopy\\CENTRIST_feature.npy\")\n",
    "X16=np.load(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\stop_sign\\CENTRIST_feature.npy\")\n",
    "X17=np.load(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\windsor_chair\\CENTRIST_feature.npy\")\n",
    "X1=np.row_stack((X11,X12,X13,X14,X15,X16,X17))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#X2(HOG FEATURE):\n",
    "X21=np.load(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\dollar_bill\\\\HOG_feature.npy\")\n",
    "X22=np.load(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\Faces\\\\HOG_feature.npy\")\n",
    "X23=np.load(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\garfield\\\\HOG_feature.npy\")\n",
    "X24=np.load(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\Motorbikes\\HOG_feature.npy\")\n",
    "X25=np.load(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\snoopy\\HOG_feature.npy\")\n",
    "X26=np.load(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\stop_sign\\HOG_feature.npy\")\n",
    "X27=np.load(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\windsor_chair\\HOG_feature.npy\")\n",
    "X2=np.row_stack((X21,X22,X23,X24,X25,X26,X27))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#X3(LBP FEATURE):\n",
    "X31=np.load(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\dollar_bill\\\\LBP_feature.npy\")\n",
    "X32=np.load(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\Faces\\\\LBP_feature.npy\")\n",
    "X33=np.load(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\garfield\\\\LBP_feature.npy\")\n",
    "X34=np.load(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\Motorbikes\\LBP_feature.npy\")\n",
    "X35=np.load(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\snoopy\\LBP_feature.npy\")\n",
    "X36=np.load(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\stop_sign\\LBP_feature.npy\")\n",
    "X37=np.load(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\windsor_chair\\LBP_feature.npy\")\n",
    "X3=np.row_stack((X31,X32,X33,X34,X35,X36,X37))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#X4(GIST FEATURE):\n",
    "X41=sio.loadmat(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\dollar_bill\\\\GIST_feature.mat\")[\"GIST_feature\"]\n",
    "X42=sio.loadmat(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\Faces\\\\GIST_feature.mat\")['GIST_feature']\n",
    "X43=sio.loadmat(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\garfield\\\\GIST_feature.mat\")['GIST_feature']\n",
    "X44=sio.loadmat(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\Motorbikes\\\\GIST_feature.mat\")[\"GIST_feature\"]\n",
    "X45=sio.loadmat(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\snoopy\\\\GIST_feature.mat\")['GIST_feature']\n",
    "X46=sio.loadmat(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\stop_sign\\\\GIST_feature.mat\")['GIST_feature']\n",
    "X47=sio.loadmat(\"F:\\lab\\code\\MLAN\\Caltech101(part)\\windsor_chair\\\\GIST_feature.mat\")['GIST_feature']\n",
    "X4=np.row_stack((X41,X42,X43,X44,X45,X46,X47))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "sio.savemat('input.mat',{'X1':X1,'X2':X2,'X3':X3,'X4':X4})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "X=[X1,X2,X3,X4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def EProjSimplex_new(v):\n",
    "    k=1\n",
    "    ft=1\n",
    "    m=v.shape[0]\n",
    "    v0=v-np.mean(v)+1/m\n",
    "    vmin=v0.min()\n",
    "    x=np.zeros((m))\n",
    "    if vmin<0:\n",
    "        f=1\n",
    "        lambda_m=0\n",
    "        while abs(f)>10**(-10):\n",
    "            v1=v0-lambda_m\n",
    "            posidx= (v1>0)\n",
    "            npos=np.sum(posidx)\n",
    "            g=-npos\n",
    "            f=np.sum(v1[posidx])-k\n",
    "            lambda_m=lambda_m-f/(g+0.1-0.1)\n",
    "            ft+=1\n",
    "            if ft>100:\n",
    "                for i in range(m):\n",
    "                    if v1[i]>0:\n",
    "                        x[i]=v1[i]\n",
    "                break\n",
    "        for i in range(m):\n",
    "            if v1[i]>0:\n",
    "                x[i]=v1[i]\n",
    "    else:\n",
    "        x=v0\n",
    "        \n",
    "    return x\n",
    "            \n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cluster(X,c,num_of_adaptive_neighbours=9):\n",
    "    \n",
    "    #Preprocess:\n",
    "    \n",
    "    n,_=X[0].shape\n",
    "    v=len(X)#v:number of different views\n",
    "    k_=num_of_adaptive_neighbours\n",
    "    lambda_=random.randint(1,30)\n",
    "    W=np.ones((v),dtype=np.float32)/v#W records the weight of each view\n",
    "    for i in range(v):\n",
    "        X[i]=pre.scale(X[i],axis=1)#You need to import sklearn.preprocessing as pre before using the standardlization module.\n",
    "    print(\"Preprocession finished...\")\n",
    "    \n",
    "    #-------------------------------------------------------\n",
    "    \n",
    "    \n",
    "    #Initialize:\n",
    "    \n",
    "    #initialize weighted_distX:\n",
    "    SUM=np.zeros((n,n),dtype=np.float32)\n",
    "    distX_initial=np.zeros((n,n),dtype=np.float32)\n",
    "    for k in range(v):\n",
    "        for i in range(n):\n",
    "            for j in range(n):\n",
    "                distX_initial[i][j]=la.norm(X[k][i]-X[k][j])**2\n",
    "        SUM=SUM+distX_initial\n",
    "    distX=1/v*SUM#distX denotes the distance matrix\n",
    "    \n",
    "    #initialize S:\n",
    "    S=np.zeros((n,n),dtype=np.float32)\n",
    "    rr=np.zeros((n,1),dtype=np.float32)\n",
    "    ascending_order_index=np.zeros((n,n),dtype=np.integer)\n",
    "    distXs=np.zeros((n,n),dtype=np.float32)\n",
    "    for i in range(n):\n",
    "        ascending_order_index[i]=np.argsort(distX[i])\n",
    "        distXs[i]=distX[i][ascending_order_index[i]]\n",
    "    for i in range(n):\n",
    "        di=distXs[i][1:k_+2]\n",
    "        rr[i]=0.5*(k_*di[k_]-np.sum(di[:k_]))\n",
    "        id=ascending_order_index[i][1:k_+2]\n",
    "        S[i][id]=(di[k_]-di)/(k_*di[k_]-np.sum(di[:k_])+10**(-10))\n",
    "    alpha=np.mean(rr)\n",
    "    \n",
    "    #initialize F:\n",
    "    S=(S+S.T)/2\n",
    "    D=np.diag(np.sum(S,axis=1))\n",
    "    L=D-S\n",
    "    eigenvalue,eigenvector=la.eig(L)\n",
    "    sorted_indices = np.argsort(eigenvalue)\n",
    "    F=eigenvector[:,sorted_indices[:c]]#F denotes the indicator matrix\n",
    "    \n",
    "    #Compute the objective value:\n",
    "    result=0\n",
    "    for k in range(v):\n",
    "        a=0\n",
    "        for i in range(n):\n",
    "            for j in range(n):\n",
    "                a+=S[i][j]*(la.norm(X[k][i]-X[k][j]))**2\n",
    "        result+=a**0.5\n",
    "    result+=alpha*(la.norm(S))**2\n",
    "    print('After initialization, the objective value is %f.'%result)\n",
    "    \n",
    "    #-------------------------------------------------------\n",
    "    \n",
    "    \n",
    "    #Update:\n",
    "    for iter_num in range(30):\n",
    "        \n",
    "        #update W:\n",
    "        for k in range(v):\n",
    "            tmp=0\n",
    "            for i in range(n):\n",
    "                for j in range(n):\n",
    "                    tmp+=S[i][j]*((la.norm(X[k][i]-X[k][j]))**2)\n",
    "            W[k]=0.5/(tmp**0.5)\n",
    "        \n",
    "        #update S:\n",
    "        S=np.zeros((n,n),dtype=np.float32)\n",
    "        ad=np.zeros((n,k_),dtype=np.float32)\n",
    "        Dist=np.zeros((n,n),dtype=np.float32)\n",
    "        for i in range(n):\n",
    "            for j in range(n):\n",
    "                Dist[i][j]+=lambda_*la.norm(F[i,:]-F[j,:])**2 \n",
    "                for k in range(v):\n",
    "                    Dist[i][j]+=W[k]*(la.norm(X[k][i,:]-X[k][j,:]))**2\n",
    "        for i in range(n):\n",
    "            idxa0=ascending_order_index[i][1:k_+1]\n",
    "            S[i][idxa0]=EProjSimplex_new(-Dist[i][idxa0]/(2*alpha))\n",
    "        \n",
    "        #update F:\n",
    "        S=(S+S.T)/2\n",
    "        D=np.diag(np.sum(S,axis=1))\n",
    "        L=D-S\n",
    "        F_old=F\n",
    "        eigenvalue,eigenvector=la.eig(L)\n",
    "        sorted_indices = np.argsort(eigenvalue)\n",
    "        F=eigenvector[:,sorted_indices[:c]]\n",
    "        \n",
    "        #update lambda_:\n",
    "        thre=10**(-10)\n",
    "        fn1=np.sum(eigenvalue[sorted_indices[:c]])\n",
    "        fn2=np.sum(eigenvalue[sorted_indices[:c+1]])\n",
    "        if fn1>thre:\n",
    "            lambda_=2*lambda_\n",
    "        elif fn2<thre:\n",
    "            lambda_=lambda_/2.0\n",
    "            F=F_old\n",
    "        else:\n",
    "            break\n",
    "            \n",
    "        #Compute the objective value:\n",
    "        result=0\n",
    "        for k in range(v):\n",
    "            a=0\n",
    "            for i in range(n):\n",
    "                for j in range(n):\n",
    "                    a+=S[i][j]*(la.norm(X[k][i]-X[k][j]))**2\n",
    "            result+=a**0.5\n",
    "        result+=alpha*(la.norm(S))**2\n",
    "        print('After iteration %i,'%(iter_num+1),end='')\n",
    "        print(' the objective value is %f.'%result)\n",
    "        \n",
    "    return S"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Preprocession finished...\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-10-6d838b0ee6d2>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mS\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mcluster\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mX\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m7\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;32m<ipython-input-9-a1758a8db058>\u001b[0m in \u001b[0;36mcluster\u001b[1;34m(X, c, num_of_adaptive_neighbours)\u001b[0m\n\u001b[0;32m     23\u001b[0m         \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mn\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     24\u001b[0m             \u001b[1;32mfor\u001b[0m \u001b[0mj\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mn\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 25\u001b[1;33m                 \u001b[0mdistX_initial\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mla\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnorm\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mX\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mX\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m**\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     26\u001b[0m         \u001b[0mSUM\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mSUM\u001b[0m\u001b[1;33m+\u001b[0m\u001b[0mdistX_initial\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     27\u001b[0m     \u001b[0mdistX\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mv\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mSUM\u001b[0m\u001b[1;31m#distX denotes the distance matrix\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mG:\\Anaconda\\setup\\lib\\site-packages\\numpy\\linalg\\linalg.py\u001b[0m in \u001b[0;36mnorm\u001b[1;34m(x, ord, axis, keepdims)\u001b[0m\n\u001b[0;32m   2254\u001b[0m                 \u001b[0msqnorm\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mreal\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mreal\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mimag\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mimag\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   2255\u001b[0m             \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 2256\u001b[1;33m                 \u001b[0msqnorm\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   2257\u001b[0m             \u001b[0mret\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0msqrt\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msqnorm\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   2258\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "S=cluster(X,7)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
